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PREFACE 

This  technical  report  is  based  on  work  performed  in  support  of  Project  20680001  at  the 
Air  Force  Armament  Laboratory,  Armament  Development  and  Test  Center,  Eglin  Air  Force 
Base,  Florida  32542.  The  work  was  performed  from  September  to  December  1974  by  the 
Systems  Analysis  and  Simulation  Branch,  Guided  Weapons  Division. 

This  technical  report  has  been  reviewed  and  is  appproved  for  publication. 
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SECTION  I 


INTRODUCTION 

The  objective  of  this  effort  is  to  develop  a package  of  programs  tor  the  Hewlett  Packard 
(HP)  9830  calculator  which  will  aid  the  engin^'er  in  the  process  of  analyzing  the  stability  and 
control  of  complex  systems  - in  particular,  guided  weapon  systems.  The  programs  presented 
here  will  aid  the  engineer  in  performing  the  linear  analysis  associated  with  continuous  systems, 
and  they  will  provide  a means  for  debugging  and  checking  more  detailed,  non-linear  simulations 
of  the  system. 

Performing  linear  analysis  on  a system  is  typically  a two-step  process.  The  first  step 
involves  the  reduction  of  a complex  block  diagram,  which  represents  the  system,  into  a simple, 
equivalent  block  diagram.  The  loop  solver  program  (LSP)  performs  this  block  diagram  reduction. 
This  program  will  accept  a block  diagram  which  contains  up  to  seven  individual  forward  transfer 
functions  and  up  to  four  feedback  transfer  functions  and  will  reduce  these  to  a single,  equiva 
lent  forward  transfer  function  and  the  associated  feedback  transfer  function. 

The  second  step  of  the  linear  analysis  involves  using  these  equivalent  transfer  functions  to 
generate  a root  locus  plot,  a Bode  stability  plot,  a Bode  response  plot,  a time  response  plot, 
or  a variety  of  other  plots  which  will  yield  information  about  the  gain  and  phase  margins,  the 
loop  frequencies  and  the  damping,  and  the  loop  response  of  the  systern.  The  stability  analysis 
program  (SAP)  will  generate  the  root  locus  plot,  the  Bode  stability  and  response  plots,  and  the 
time  response  plot  of  the  system. 

There  are  two  versions  of  the  stability  analysis  executive  program.  One  version  is  designed 
to  take  the  equivalent  transfer  functions  directly  from  the  LSP  by  means  of  a data  tape.  This 
version  is  referred  to  as  Stability  Analysis  Program-Tape  Input  (SAP-TI).  The  second  is  designed 
to  allow  the  operator  to  enter  a single  forward  transfer  function  and  a single  feedback  transfer 
function  by  means  of  the  keyboard,  without  having  to  resort  to  the  LSP.  This  second  version 
is  referred  to  as  the  Stability  Analysis  Program-Keyboard  Input  (SAP-KI). 

Due  to  the  size  of  the  second  version  of  the  SAP,  the  entire  program  cannot  reside  in 
memory  at  one  time.  Therefore,  each  SAP  consists  of  an  executive  routine,  which  resides  in 
memory,  and  three  subroutines,  only  one  of  which  resides  in  memory  at  any  one  time.  The 
basic  interaction  of  these  programs  is  shown  in  Figure  1. 

This  report  describes  the  loop  solver  program  and  the  two  versions  of  the  stability  analysis 
program,  along  with  the  three  subroutines  and  their  use.  A discussion  of  block  diagram  reduc 
tion  technique  and  stability  analysis  techniques  is  contained  in  Reference  1,  and  the  HP  9830 
BASIC  language  is  presented  in  Reference  2. 


References: 

1.  D'Azzo  and  Houpis,  Feedback  Control  System  Analysis  and  Synthesis.  McGraw-Hill,  N.  Y., 
1966. 

2.  Hewlett-Packard  9830A  Calculator  Operating  and  Programming  Manual.  Hewlett-Packard 
Company,  Loveland,  Colorado,  1973. 


Figure  1.  Program  Interaction  Chart 
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PROGRAM  REQUIREMENTS  AND  RESTRICTIONS 

1.  REQUIREMENTS 

The  programs  described  in  this  report  require  a HP  9830  calculator  with  the  following 
options: 

' HP  9862A  Calculator  Plotter 

Plotter  ROM 
String  Variable  ROM 
Matrix  ROM 
8K  words  of  memory 

2.  RESTRICTIONS 

The  restrictions  enumerated  below  apply,  in  general,  to  all  of  the  programs  and  their  opera 
tion.  In  those  instances  where  the  restriction  applies  to  a single  program,  it  is  so  noted. 

a.  The  system  block  diagram  must  be  reduced  to  the  form  shown  in  Figure  2 prior  to  being 
entered  into  the  loop  solver  program. 

b.  The  data  storage  file  length  must  be  at  least  2000  words. 

c.  The  order  of  the  total  system  being  processed  must  not  exceed  10th.  (Note:  The  order 
is  determined  by  adding  the  order  of  the  denominator  of  each  forward  and  feedback  transfer 
function.) 

d.  The  programs  assume  negative  feedback.  If  the  feedback  is  positive,  a negative  sign  must 
be  entered  as  a part  of  the  feedback  transfer  function. 

e.  All  polynomials  are  entered  from  the  lowest  order  coefficient  to  the  highest  order 
coefficient.  Constants  are  entered  as  polynomials  of  order  zero. 

f.  In  the  loop  solver  program,  all  transfer  functions  are  initialized  to  zero.  However,  if 

it  is  necessary  to  zero  out  a non-zero  transfer  function,  enter  a zero  in  the  numerator  and  a one 
in  the  denominator.  The  program  does  not  recognize  0/0!  This  restriction  also  applies  to  enter 
ing  a feedback  transfer  function  of  zero  in  the  keyboard  input  version  of  the  SAP-KI. 

I g.  When  the  calculator  asks  a question  which  requires  a letter  response,  e.g.,  "Y"  or  "N" 

for  "yes"  or  "no",  the  appropriate  letter  must  be  entered.  Any  other  input,  including  the 
command  "STOP",  will  result  in  an  error,  i.e.,  ERROR  77.  This  error  is  a recoverable  error, 

» and  execution  will  continue  when  the  appropriate  letter  is  entered. 
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h.  The  case  identifier  must  not  exceed  ten  alphanumeric  characters. 


i.  When  using  either  of  the  two  versions  of  the  stability  analysis  programs  and  when  the 
execution  is  under  the  control  of  one  of  the  three  subroutines,  halting  the  program  and  then 
continuing  it  will  result  in  an  error.  This  is  due  to  the  fact  that  the  RETURN  flag  is  lost 
when  the  program  is  halted.  The  way  to  recover  from  this  error  is  to  stop  the  execution  and 
continue  execution  at  some  step  in  tfie  main  stability  analysis  program. 
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LOOP  SOLVER  PROGRAM  (LSP) 

1.  PROGRAM  DESCRIPTION 

The  loop  solver  program  (LSP)  has  been  designed  to  accept  a complex  block  diagram  with 
up  to  seven  forward  transfer  functions  and  four  feedback  transfer  functions,  and  to  reduce  them 
to  an  equivalent  single  forward  transfer  function  and  an  associated  equivalent  feedback  transfer 
function.  The  input  data  and  the  calculated  equivalent  transfer  functions  are  stored  on  magnetic 
tape  and  are  printed  out.  Data  tape  storage  allows  the  transfer  of  the  resulting  polynomials  to 
the  tape  input  version  of  the  SAP-TI. 

The  LSP  is  set  up  to  process  a block  diagram  of  the  form  shown  in  Figure  2.  Each  of  the 
transfer  functions  may  consist  of  several  polynomials  which  will  be  multiplied  together  as  they 
are  entered  into  the  program.  All  of  the  individual  transfer  functions  are  initialized  to  zero, 
so  it  is  only  necessary  to  enter  non-zero  functions.  The  total  system  must  not  exceed  10th 
order. 


The  LSP  will  generate  two  equivalent  transfer  functions  for  each  of  the  nodes  shown  in 
Figure  2.  These  two  transfer  functions  are  designated  as  (1)  an  internal  loop  equivalent  and 
(2)  an  isolated  loop  equivalent.  This  terminology  is  not  standard;  thus,  these  terms  are  defined 
as  follows: 


Internal  Loop  Equivalent  - The  internal  loop  equivalent  is  the  transfer  function  that  is 
generated  when  all  transfer  functions  between,  or  internal  to,  the  designated  node,  e.g.,  2,  and 
its  mirror  node,  e.g.,  2',  are  reduced  to  a single  transfer  function. 

Isolated  Loop  Equivalent  - The  isolated  loop  equivalent  is  the  transfer  function  that  is 
generated  when  all  transfer  functions,  both  internal  to  and  external  to  the  designated  node  and 
its  mirror  node,  are  reduced  to  a single  transfer  function. 


The  two  types  of  equivalents  are  illustrated  for  node  2 in  Figure  3.  As  can  be  seen  in 
Figure  3,  the  isolated  loop  equivalent  is  a combination  of  the  internal  loop  equivalent  plus 
an  equivalent  feedback,  which  is  due  to  the  external  transfer  functions.  In  general. 


X “ ^x-1 


1 + H^GI^.i 


G] 


x+1 


where 

G]  ^ = the  internal  loop  equivalent  transfer  function  for  node  x 
G^  = the  forward  transfer  function  numbered  x 

= the  feedback  transfer  function  that  is  fed  into  node  x 
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The  effect  of  the  transfer  functions  external  to  node  x can  be  considered  as  a single  feedback 
and,  in  general,  is 


Hl^=  + G4+)(G4  X 


where 


H]  ^ = the  equivalent  feedback  that  is  external  to  node  x,  i.e.,  the  feedback  directly 
into  node  x is  not  included 

Having  the  equivalent  forward  transfer  function,  G]j^,  and  the  equivalent  feedback  transfer 
function,  H]  for  node  x,  the  isolated  equivalent  transfer  function  is  then  defined  as 


G], 


■ 1 + GU  HI, 


where 

G]  * = the  isolated  loop  equivalent  transfer  function  for  node  x 

It  should  be  noted  that,  at  node  4,  the  internal  and  isolated  equivalents  are  equal  since  there 
are  no  loops  external  to  node  4,  and  hence,  the  equivalent  external  feedback,  i.e.,  H]4,  is 
zero. 

Normally,  the  internal  loop  equivalents  would  be  used  when  designing  a servo  system. 

The  isolated  loop  equivalents  are  useful  when  a system  has  been  designed  and  one  wants  to 
examine  the  sensitivity  of  a single  loop  in  the  presence  of  the  rest  of  the  system.  The  isolated 
loop  equivalents  are  also  useful  when  using  the  results  of  the  linear  analysis  to  verify  full-up 
system  simulation. 

However,  caution  must  be  used  in  interpreting  the  analysis  results  when  using  the  isolated 
Inrjp  equivalent  forward  transfer  functions.  This  is  due  to  the  way  in  which  the  equivalent 
foiward  transfer  functions  are  derived.  The  analysis  subroutines  associated  with  SAPs  assume 
that  the  closed-loop  system  transfer  function  is  of  the  form 


C-  = KG 
R 1 -t-  KGH 

The  gain  that  is  entered  as  a part  of  the  analysis  subroutines  is  the  K in  the  above  equation  and 
is  associated  with  the  forward  transfer  function,  G,  or  in  this  case,  the  equivalent  forward  trans- 
fer function,  Glx  or  G]  *.  The  fact  that  the  gain  is  associated  with  the  isolated  loop  equivalent 
forward  transfer  function  means  that  the  gain  is  distributed  in  the  individual  forward  transfer 
functions  and  in  the  internal  feedback  transfer  functions,  i.e.,  the  gain  cannot  be  associated 
with  any  one  single  transfer  function  in  the  original  system.  This  limits  the  utility  of  the 
isolated  loop  equivalent  in  those  cases  where  the  closed-loop  transfer  function  is  used,  e.g., 
the  closed  loop  Bode  plot  and  the  time  response  programs. 
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2.  OPERATING  PROCEDURE 


a.  Load  the  LSP  into  the  HP  9830. 

b.  Insert  a tape  with  a 2000  word  data  file  in  the  calculator. 

c.  Press  RUN/EXECUTE.  The  display  will  request  the  number  of  the  forward  transfer 
function  that  is  to  be  entered.  An  entry  of  zero  will  cause  the  program  to  continue  to  the 
feedback  transfer  functions  (step  f).  If  a non-zero  entry  is  made,  the  printer  will  print  "ENTER 
NUMERATOR  DATA  FOR  G#  (entry)"  and  the  display  will  request  the  order  of  the  first 
polynomial  in  the  numerator. 

d.  After  the  order  is  entered,  the  display  will  request  the  coefficients  starting  with  tfie 
lowest  order  coefficient.  A zero  order  polynomial  with  the  constant  equal  to  zero  will  termi- 
nate the  input  of  the  numerator. 

e.  The  printer  will  print  "ENTER  DENOMINATOR  DATA  FOR  G#  (entry)",  and  the 
display  will  ask  for  the  order  of  the  first  polynomial  in  the  denominator.  A zero  order  poly 
nomial  with  the  constant  equal  to  zero  will  cause  the  program  to  return  to  step  c. 

f.  The  display  will  request  the  number  of  the  feedback  transfer  function  that  is  to  be 
entered.  An  entry  of  zero  will  cause  the  program  to  continue  to  step  g.  A non-zero  entry 
allows  the  entry  of  data,  as  in  steps  d and  e. 


g.  The  forward  and  feedback  transfer  function  polynomials  will  be  printed  out  in  matrix 
format.  The  number  of  the  column  in  each  matrix  corresponds  to  the  number  of  the  transfer 
function,  as  defined  in  Figure  2.  The  first  row  of  the  matrix  is  the  zero  order  coefficient  and 
so  on  up  to  the  10th  order  coefficient  which  is  in  row  11.  The  12th  row  contains  the  order 
of  the  polynomial  which  is  contained  in  that  column. 

h.  The  calculator  will  then  calculate  the  equivalent  loop  transfer  functions.  This  is  time 
consuming  , and  the  calculator  will  appear  to  be  totally  inactive  during  this  time. 

i.  The  calculated  internal  and  isolated  loop  equivalent  transfer  functions  will  be  printed 
in  the  same  format  as  in  step  g.  The  number  of  the  columns  corresponds  to  the  number  of 
the  node,  as  defined  in  Figure  2. 

j.  The  display  will  request  the  number  of  the  data  file  in  which  the  calculated  data  is  to 
be  stored.  After  the  data  is  stored,  the  program  will  return  to  step  c.  This  allows  one  to 
easily  store  a new  set  of  equivalent  loop  data  for  a system  which  has  only  minor  differences 
from  the  original  one.  Since  the  previous  data  for  all  the  forward  and  feedback  transfer  functions 
are  still  in  memory,  it  is  only  necessary  to  enter  those  which  are  different  in  the  new  system. 

If  this  option  is  not  desired,  merely  press  STOP. 

3.  EXAMPLE  PROBLEM 

The  following  system  block  diagram  (Figure  4)  can  be  put  into  the  generalized  format 
(Figure  2)  in  several  ways.  Two  ways  will  be  illustrated  here: 
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liMEF  DEHliMINFiTOP  1h;Th  FOR  Utt  1 
FHTFR  HUNERR  f OR  iiHTR  F i.iR  Fi# 
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FirrEF  HIJMEFRTGR  PRTfi  FOR  G#  3 
ENTER  liENOFUNRiOR  URTh  FOR  G# 

E/I!ER  NUMERRTOF  URTH  FOR  G#  4 
EIFTER  DEHOMiyR'GF:  BRTR  FOR  G#  4 
FMTER  NLlTlEPRTW  DRTP  FOR  G#  5 
ENTER  CEMOFiTHfiTOR  riRTR  FOR  G#  G 
I.MTER  HUilERrt  TOR  PR  FR  FOR  G#  R. 

ENTER  PEHOMIHRTOR  LiRTR  FOR  G#  6 
ENTER  NUMERATOR  IiRTFl  FOR  G#  7 
ENTEF’  riEHOl'nNRTOR  DRTR  FOR  Gtt  7 
INTER  NUMERRiOR  DRTR  FOR  H#  3 
EN-;tT?«''iSN0"lIKV70P-  FQP-HN  3 

Figure  5.  Input  Sequence  for  LSP  (Example  1) 
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Figure  6.  LSP  Output  (Example  1) 
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Figure  6.  LSP  Output  (Example  1)  (Concluded) 


Example  2: 


Let  G4  = (s  t.  1L.(0-5)  (5QQ) 

(s  + 10)  (s)  (s  + 10)  (s2),r 

and 


Hi  = (s+  1)/1 

The  sequence  in  which  the  data  was  entered  is  shown  in  Figure  7,  and  the  results  are  shown 
in  Figure  8. 

It  can  be  seen  that  the  data  for  node  1 in  example  2 is  exactly  equal  to  the  data  for 
node  3 in  example  1. 

4.  FLOW  CHART 

The  flow  chart  for  the  LSP  is  shown  in  Figure  9.  The  numbers  in  parentheses  on  the  right 
and  left  side  of  the  chart  are  the  line  numbers  that  should  be  used  to  initiate  that  particular 
function  or  process  in  the  flow  diagram.  The  complete  listing  for  the  program  is  contained 
in  Appendix  A. 
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Figure  7.  Input  Sequence  for  LSP  (Example  2 
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SECTION  IV 


STABILITY  ANALYSIS  PROGRAM  TAPE  INPUT 
(SAP-TI) 


1.  PROGRAM  DESCRIPTION 

The  SAP-TI  has  been  designed  to  retrieve  designated  transfer  functions  from  the  data  tafie 
generated  by  the  LSP  and  to  control  the  execution  of  the  three  analysis  subroutines,  i.e.,  the 
Root  Locus  (RL),  the  Bode  (B),  and  the  Time  Response  (TR)  subroutines. 

The  program  will  load  common  with  data  from  the  user  specified  data  file  and  then 
retrieve  the  equivalent  forward  and  feedback  transfer  functions  for  the  node  and  loop  type 
which  are  specified  by  the  user.  The  program  then  constructs  three  polynomials  which  are  used 
in  the  three  subroutines.  The  three  polynomials  are  as  follows; 

NG«DH  - (Numerator  of  G) ‘(Denominator  of  H) 

NG*NH  = (Numerator  of  G)*(Numerator  of  H) 

DG*DH  = (Denominator  of  G)*(Denominator  of  H) 

Both  the  characteristic  equation  and  the  closed  loop  transfer  function  can  be  obtained  from  these 
three  polynomials  since 

C = KG  = KING)  (DH> 

R 1 -t  KGH  (DGKDH)  + K(NG)(NH) 

The  program  tnt..i  loads  the  appropriate  subroutine  and  transfers  control  to  the  subroutine. 

2 OPERATING  PROCEDURE 

a.  Load  the  SAP-TI  into  the  HP  9830. 

b.  Insert  the  data  tape  in  the  calculator. 

c.  Press  RUN/EXECUTE.  The  display  will  request  the  number  of  the  data  file.  The  data 
will  be  loaded  into  common  from  the  designated  data  file. 

d.  The  display  will  request  the  loop  type,  i.e.,  internal  or  isolated.  The  definitions  of 
the  internal  and  isolated  loops  are  the  same  as  in  the  LSP  description. 

e.  The  display  will  request  the  node  number  which  is  defined  in  Figure  2. 

f.  The  printer  will  print  the  statement  "ROOT  LOCUS  (L),  BODE  (B),  RESPONSE  (R), 

OR  STOP  (S)?".  There  will  be  a short  time  delay  between  steps  e and  f because  the  calculator 
is  generating  the  three  polynomials  referred  to  in  the  program  description.  Before  entering 

the  proper  alphanumeric  character  in  response  to  the  printed  question,  be  sure  that  the  program 
tape  is  replaced  in  the  calculator. 


19 


1 


g.  The  calculator  will  load  the  selected  subroutine  into  memory.  There  are  internal  flags 
set  in  the  program  so  that  this  step  will  be  skipped  if  the  selected  subroutine  is  already  in 
memory. 

t.  h.  At  this  point,  the  selected  subroutine  will  be  executed.  The  operating  procedures  for 

the  subroutines  are  contained  in  the  following  sections  of  this  report. 

i.  After  the  subroutine  is  executed,  the  display  will  ask  if  a new  transfer  function  is 
desired.  If  the  response  is  negative,  the  program  will  return  to  step  f. 

j.  If  a new  transfer  function  is  requested,  the  display  will  ask  whether  the  new  transfer 
function  is  in  the  existing  data  file  or  in  a new  data  file.  If  the  existing  file  is  indicated,  the 
program  will  return  to  step  d. 

k.  If  the  new  transfer  function  is  in  a new  data  file,  put  the  data  tape  which  contains 
the  appropriate  data  file  in  the  calculator.  The  display  will  request  rhe  number  of  the  data 
file,  the  data  will  be  loaded  into  common,  and  the  program  will  return  to  step  d. 

3.  FLOW  CHART 

^ The  flow  chart  for  the  SAP-TI  program  is  shown  in  Figure  10.  The  numbers  in  parentheses 

on  the  right  and  left  sides  of  the  chart  are  the  line  numbers  that  should  be  used  to  continue 
program  execution  at  that  particular  function  or  process  on  the  flow  diagram.  The  complete 
listing  for  the  program  is  contained  in  Appendix  B. 

Two  flags  are  used  to  keep  track  of  which  of  the  three  subroutines  are  in  memory.  Two 
flags,  one  in  common  and  one  in  memory,  are  used  because  loading  a new  data  file  will  destroy 
the  hag  that  is  in  common  and  the  loading  of  a new  subroutine  will  destroy  the  flag  that  is 

1 not  in  common.  The  use  of  the  two  flags  insures  that  there  is  always  one  flag  to  indicate  which 

E subroutine  is  in  memory. 

1 

i 


1 
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SECTION  V 


[ 

I 


STABILITY  ANALYSIS  PROGRAM-KEYBOARD  INPUT 
(SAP  Kl) 

1.  PROGRAM  DESCRIPTION 

The  SAP-KI  is  a combination  of  the  LSP  and  the  SAP-TI  programs.  The  program  allows 
the  user  to  enter  the  forward  and  feedback  functions  directly  into  the  stability  analysis  program 
without  having  to  use  the  LSP  and  data  tape.  This  program  is  useful  when  the  system  is  simple 
enough  that  the  block  diagram  reduction  can  be  done  manually. 

2.  OPERATING  PROCEDURE 

a.  Load  the  SAP-KI  into  the  HP  9830. 

b.  Press  RUN/EXECUTE.  The  printer  will  print  "ENTER  NUMERATOR  OF  G",  and 
the  display  will  request  the  order  of  the  first  polynomial  in  the  numerator. 

c.  After  the  order  is  entered,  the  display  will  request  the  coefficients  starting  with  the 
lowest  order  coefficient.  A zero  order  polynomial  with  the  numerator. 

d.  The  printer  will  print  "ENTER  DENOMINATOR  OF  G",  and  the  display  will  request 
the  order  of  the  polynomial  as  in  step  b and  the  coefficients  are  entered  as  in  step  c. 

e.  The  numerator  of  H is  entered  as  in  steps  b and  c. 

f.  The  denominator  of  H is  entered  as  in  steps  b and  c. 

g.  The  printer  will  print  the  statement  "ROOT  LOCUS  (L),  BODE  (B),  RESPONSE  (R), 
OR  STOP  (S)?".  There  will  be  a short  time  delay  between  steps  f and  g because  the  calcu 
lator  is  generating  the  three  polynomials  referred  to  in  the  SAP-TI  Program  Description  in 
Section  III. 

h.  The  calculator  will  load  the  selected  subroutine  into  memory.  There  are  internal  flags 
set  in  the  program  so  that  this  step  will  be  skipped  if  the  selected  subroutine  is  already  in 
memory. 

i.  At  this  point,  the  selected  subroutine  will  be  executed.  The  operating  procedures  for 
the  subroutines  are  contained  in  the  following  sections  of  this  report. 

j.  After  the  subroutine  is  executed,  the  display  will  ask  if  a new  transfer  function  is 
desired.  If  the  response  is  negative,  the  program  will  return  to  step  g. 

j.  If  a new  transfer  function  is  requested,  the  program  will  return  to  step  b. 

3.  EXAMPLE  PROBLEM 

Figure  11  shows  the  keyboard  input  sequence  for  the  system  shown  in  Figure  4.  For  this 
example  the  PRINT  ALL  option  on  the  calculator  was  used  to  show  the  information  that  is 
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Figure  11.  Input  Sequence  for  SAP-KI 


normally  on  the  display.  In  normal  operation  only  the  four  statements  that  are  designated  hy 
the  arrow  would  appear  on  the  printer.  The  questions  will  normally  appear  on  the  display, 
and  the  number  following  the  question  mark  is  the  number  that  the  user  inputs  in  response  to 
that  question. 

4.  FLOW  CHART 

The  flow  chart  for  the  SAP-KI  program  is  shown  in  Figure  12.  The  numbers  in  parantheses 
on  the  right  and  left  sides  of  the  chart  are  the  line  numbers  that  should  be  used  to  continue 
the  program  execution  at  that  particular  function  or  process  on  the  flow  diagram.  The  complete 
listing  for  the  program  is  contained  in  Appendix  C.  F4  is  a flag  to  indicate  what  subroutine 
is  in  memory. 


SECTION  VI 


ROOT  LOCUS  SUBROUTINE 

1.  SUBROUTINE  DESCRIPTION 

The  Root  Locus  (RL)  Subroutine  is  used  in  conjunction  with  either  the  SAP-TI  executive 
program  or  the  SAP-KI  executive  program.  This  subroutine  finds  the  loci  of  the  roots  of  the 
characteristic  equation  associated  with  the  system  transfer  function.  The  subroutine  also  pro 
vides  a listing  of  the  roots  and  a plot  of  the  loci  if  the  user  desires. 

The  subroutine  forms  the  characteristic  equation  with  two  of  the  three  polynomials  referenced 
in  Section  IV.  The  resulting  form  of  the  characteristic  equation  is 

(DG)(DH)  + K (NG)(NH)  = 0 
where  K is  the  variable  gain  in  the  system. 

The  zeros  are  the  roots  of  the  term  (NG)(NH),  and  the  poles  are  the  roots  of  the  term  (DG) 
(DH).  The  closed  loop  poles  are  the  roots  of  the  characteristic  equation  with  a specified  value 
of  the  gain,  K. 

2.  OPERATING  PROCEDURE 

a.  The  display  will  request  a CASE  IDENTIFIER.  Type  in  up  to  ten  alphanumeric  characters 
that  you  wish  to  use  to  identify  the  following  plots  and  print-outs.  After  the  identifying  charac- 
ters, press  either  EXECUTE  or  STOP  to  continue  program  execution. 

b.  The  case  identifier,  the  numerator  coefficients,  and  the  denominator  coefficients  will  be 
printed. 

c.  The  display  will  ask  whether  a plot  is  desired.  If  the  response  is  negative,  the  program 
will  continue  with  step  f. 

d.  If  a plot  was  requested,  the  display  will  request  the  maximum  values  for  the  real  and 
imaginary  axes  and  for  the  increments  on  these  two  axes.  The  sign  associated  with  the  values 
input  is  not  important  since  the  calculator  only  uses  the  magnitude  and  affixes  the  proper 
sign.  The  program  draws  a six-inch  real  axis  and  a five-inch  imaginary  axis,  so  the  plot  is  more 
readable  if  the  maximum  value  of  the  real  axis  is  a multiple  of  six  and  the  maximum  imaginary 
value  is  a multiple  of  five.  Prior  to  entering  the  imaginary  increment,  be  sure  that  the  plotter 
is  set  up  properly. 

e.  The  axes  will  be  drawn  and  labeled.  It  should  be  noted  that  the  program  only  plots 
one  quadrant  of  the  S-plane. 

f.  The  zeros  will  be  calculated  and  printed  (and  plotted).  The  first  number  printed  is 
the  real  part  of  the  root,  and  the  second  number  is  the  imaginary  part  of  the  root. 
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g.  The  poles  will  be  calculated  and  printed  (and  plotted). 


I. 


h.  The  display  will  request  the  minimum  gain,  maAimum  gain,  and  gam  incren.ent.  The 
gain  is  the  value  of  K referenced  in  the  Program  Description,  alwve, 

i.  The  roots  of  the  characteristic  equation  for  each  value  of  gain  which  was  specified  m step 

^ h will  be  calculated  and  printed  (and  plotted). 

j.  The  display  will  ask  if  a new  gain  increment  is  desired.  If  the  response  is  affirmative, 

C the  program  will  return  to  step  h.  If  the  response  is  negative  and  if  a plot  was  not  originally 

requested,  the  program  will  return  control  to  the  SAP  executive  program. 

k.  If  the  response  is  negative  and  if  a plot  was  originally  requested,  the  display  will  ask 
if  a replot  is  desired.  This  option  is  provided  in  case  the  user  wants  to  expand  some  portion 

. of  the  original  plot.  If  a plot  was  originally  requested,  the  coordinates  of  each  point  that  was 

plotted  were  also  stored  in  a storage  array.  Due  to  the  dimensions  of  the  storage  array,  only 
the  first  65  points  are  stored.  When  a replot  is  requested,  the  coordinates  of  each  previously 
stored  point  are  retrieved  from  the  storage  array,  examined  to  see  if  they  are  within  the  new 
scale,  and  if  so,  the  point  is  plotted.  If  the  response  to  the  request  for  a replot  is  negative, 
the  program  control  will  be  returned  to  the  SAP  executive  program.  If  the  response  is 
affirmative,  the  display  will  request  information  about  the  replot  axes  as  in  steps  d and  e, 
above.  Prior  to  entering  the  value  of  the  imaginary  increment,  be  sure  that  the  plotter  is  set 
up. 

1.  The  previously  calculated  roots  will  be  r^olotted;  they  will  not  be  recalculated.  The 
•.  program  will  then  return  to  step  k,  above. 

3.  EXAMPLE  PROBLEM 

The  system  in  Figure  4 is  used  in  this  example.  The  SAP-KI  executive  program  is  used. 

The  printer  output  is  shown  in  Figure  13.  The  option  for  a NEW  GAIN  INCREMENT  was 
used.  The  original  gain  increment  was  from  0.1  to  1 in  steps  of  0.1.  The  additional  gain  incre- 
ment was  from  2 to  16  in  steps  of  2.  The  results  of  this  option  can  be  seen  in  the  printout 
and  in  the  spacing  of  the  roots  on  the  plot.  The  original  plot  is  shown  in  Figure  14,  and  a 
sample  replot  is  shown  in  Figure  15.  Note  the  difference  in  the  value  of  the  axes  in  the  two 
plots.  The  replot  greatly  expands  the  movement  of  the  roots  in  the  vicinity  of  the  origin. 

4.  FLOW  CHART 

The  flow  chart  for  the  RL  subroutine  is  shown  in  Figure  16.  The  numbers  in  parentheses 
on  the  right  and  left  of  the  chart  are  the  line  numbers  that  should  be  used  to  initiate  that  parti- 
cular function  or  process  in  the  flow  diagram.  The  complete  listing  for  the  subroutine  is 
contained  in  Appendix  D. 
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Figure  13.  RL  Subroutine  Output 
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Figure  14.  RL  Subroutine  Plot 
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rigure  15.  RL  Subroutine  Replot 
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SECTION  VII 


BODE  SUBROUTINE 


1.  PROGRAM  DESCRIPTION 

The  Bode  (B)  Subroutine  is  used  in  conjunction  with  either  the  SAP  TI  executive  (program 
or  the  SAP-KI  executive  program.  The  subroutine  will  calculate  and  plot  either  the  Bode 
stability  criteria  or  the  Bode  frequency  response. 

The  Bode  stability  criteria,  which  is  referred  to  as  the  "open  loo()"  response  in  this  program, 
are  based  on  the  term  (GH)  of  the  characteristic  equation.  A plot  of  the  amplitude  and  phase 
of  the  term  (GH)  will  yield  the  gain  and  phase  margins  of  the  system,  in  this  portion  of  the 
program,  the  gain  is  taken  as  unity.  A gain  other  than  one  will  simply  translate  the  amplitude 
curve  in  the  Y direction.  The  phase  margin  is  defined  as 

PM  = 180  - j phase  I evaluated  at  that  frequency 

where  the  amplitude  is  one. 

The  gain  margin  is  defined  as 

GM  = |amplitude(  evaluated  at  that  frequency 

where  the  phase  is  180  degrees  . 

The  Bode  frequency  response  which  is  referred  to  as  the  "closed  loop"  response  in  this  pro- 
gram is  based  on  the  closed  loop  transfer  function 

C =,  KG 

R 1 -I-  K (GH) 

In  this  case,  the  response  of  the  system  is  dependent  on  the  gain  of  the  system,  and,  therefore, 
a provision  for  entering  a gain  is  incorporated  in  the  program.  As  implied  by  the  sign  in  the 
denominator  of  the  above  equation,  the  program  assumes  negative  feedback.  A system  which 
has  positive  feedback  can  be  evaluated  by  affixing  a negative  sign  to  the  feedback  transfer 
function. 

2.  OPERATING  PROCEDURE 

a.  The  display  will  request  a CASE  IDENTIFIER.  Type  in  up  to  ten  alphanumeric  charac 
ters  that  you  wish  to  use  to  identify  the  following  plots  and  print  outs.  After  the  identifying 
characters  are  entered,  press  either  EXECUTE  or  STOP  to  continue  program  execution. 

b.  The  case  identifier  will  be  printed,  and  the  display  will  ask  whether  the  open  or  closed 
loop  version  is  desired.  If  the  response  is  open  loop,  the  program  will  proceed  to  step  d. 

c.  If  the  response  is  closed  loop,  the  display  will  request  the  desired  gain. 
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d.  The  numerator  coefficients  and  the  denominator  coefficients  will  be  printed,  and  the 
display  will  request  the  frequency  range  over  which  the  response  is  to  be  evaluated  ’ The  range 
is  in  powers  of  ten.  For  example,  if  the  desired  frequency  range  is  0.1  rad/sec  to  100  rad/sec 
the  input  to  the  program  would  be  1,  2 since  C.1  = 10  ’ and  100  = 10^. 

e.  The  calculatoi  will  then  calculate  the  amplitude  and  phases  for  the  specified  frequency 
range.  This  is  time  consuming,  and  the  calculator  will  appear  to  be  totally  inactive  for  several 
minutes. 

f.  The  display  will  request  whether  a listing  of  the  phase  and  amplitude  is  desired.  A listing 
will  or  will  not  be  printed,  depending  on  the  response.  The  program  will  automatically  scale 
the  amplitude  plot  up  to  a maximum  of  120  dB.  The  phase  is  scaled  from  +180  degrees  to 

130  degrees.  The  relative  magnitude  and  the  quadrant  of  the  phase  angle  are  correct,  but  the 
absolute  magnitude  of  the  phase  angle  will  have  to  be  determined  by  the  user.  For  example 
in  the  sample  problem,  the  phase  angle  for  very  low  frequencies  should  be  approximately 
-270  degrees,  i.e.,  two  lead  terms  and  five  lag  terms.  The  program  will  only  determine  that  the 
phase  angle  due  to  the  five  lags  is  approximately  90  degrees  and  is  in  the  first  quadrant.  The 
result  of  this  is  that  the  total  phase  angle  at  low  frequencies  will  be  calculated  as  +90  degrees 
instead  of  the  correct  -270  degrees. 

g.  If  the  open  loop  response  was  originally  requested,  the  gain  and  phase  margins  will  be 
printed  before  program  control  is  returned  to  the  SAP.  This  printout  will  not  occur  if  the 
closed  loop  response  was  originally  requested. 

3.  EXAMPLE  PROBLEM 

The  system  in  Figure  4 is  used  in  this  example.  The  open  loop  response  (Bode  stability 
criteria)  printout  is  shown  in  Figure  17,  and  the  corresponding  plot  is  shown  in  Figure  18. 

For  the  closed  loop  response,  a gain  of  6.5  was  used  since  the  root  locus  plot  (Figures  14  and 
15)  indicated  that  this  gain  in  the  forward  loop  would  give  an  adequate  system  response.  The 
printout  for  hte  closed  loop  response  (Bode  frequency  response)  is  shown  in  Figure  19  and  the 
corresponding  plot  is  shown  in  Figure  20. 

4 FLOW  CHART 

The  flow  chart  for  the  Bode  subroutine  is  shown  in  Figure  21.  The  complete  program 
listing  for  the  subroutine  is  contained  in  Appendix  E. 
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Figure  19.  3 Subroutine  Output  for  Ciose'f  -Loop  Response 
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SECTION  VIII 

TIME  RESPONSE  SUBROUTINE 


1.  PROGRAM  DESCRIPTION 

The  Time  Response  (TR)  Subroutine  is  used  in  conjunction  with  either  the  SAP  TI  execu- 
tive program  or  the  SAP-KI  executive  program.  This  subroutine  will  calculate  and  plot  the 
time  response  of  a system  to  a unit  step  input.  The  time  response  is  calculated  from  inverse 
Laplace  transform.  A description  of  the  theory  used  to  calculate  the  inverse  Laplace  trans- 
form is  contained  in  Appendix  G and,  consequently,  will  not  be  discussed  here. 

The  subroutine  is  limited  in  that  the  subroutine  cannot  handle  a system  which  has  a real 
root  that  is  repeated  more  than  three  times  or  a repeated  complex  root  pair.  Also  there  is  no 
provision  for  altering  the  driving  function  for  the  system. 


The  subroutine  prints  out  information  which  will  allow  the  user  to  construct  the  time 
response  equation,  e.g.,  the  inverse  Laplace  transform,  for  the  system.  This  information  is 
labeled  "Inverse  Laplace  Transform  Information,"  and  an  example  output  is  shown  in  Figure 
22.  A FORM  1 term  is  written  as 

y (t)  = (A  -r  Bt  + Ct^je^^^ 

and  a FORM  2 term  is  written  as 

y(t)  = A cos  It  -H  B sin  lt)e^^ 

The  total  time  response  equation  is  the  sum  of  the  particular  solution  plus  the  homogeneous 
terms  listed  on  the  output.  The  particular  solution  is  always  of  the  form 


where 


Yp(t) 


aq  q^ 


q = the  order  of  the  lowest  order,  non-zero  coefficient  in  the  denominator  of  the 
system  transfer  function. 


For  example,  if  the  system  has  a pole  at  the  origin,  i.e.,  ag  = o and  a-]  4 o,  then 

YJt)  = t 

^1 

The  total  time  response  equation  for  the  example  case  shown  in  Figure  23  is 


t(t)  = 0.418bt':  - I 0.5512  cos  (1.273)t  + 0.3172  sin  (1.273)t  | e 

0.0099  ^ o.03805e  + o.007745e 


Subroutine  Output  (Example 


CH5E  EXAMPLE 
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44 


s^x^ssaasas 


Figure  23.  TR  Subroutine  Plot  (Example  1,  Ga 


The  plot  of  ttie  time  response  is  autoinatically  scaled  and  labeled.  This  is  to  insure  that 
the  plot  is  readable.  To  accomplish  this  scaling,  tfie  program  will  (1)  determine  tlie  maximum 
value  of  the  time  response  of  the  system  in  the  time  interval  specified  by  the  user,  and  12) 
determine  if  the  response  has  some  final  steady  state  value^.  If  tfie  system  does  have  a final 
steady  state  value,  tfie  program  will  scale  the  plot  as  a percentage  of  that  final  steady  state 
value,  and  the  y axis  will  be  labeled  RESPONSE  (%).  If  the  system  does  not  have  a steady  state 
value,  the  program  will  scale  the  plot  based  on  the  maximum  value  of  the  response,  and  the 
y axis  will  be  labeled  RESPONSE  (UNITS).  In  either  case,  the  scale  determined  by  the  program 
will  be  labeled  on  the  plot. 

2.  OPERATING  PROCEDURE 

a.  The  display  will  request  a CASE  IDENTIFIER.  Type  in  up  to  ten  alphanumeric 
characters  that  you  wish  to  use  to  identify  the  following  plots  and  printouts.  After  the 
identifying  characters  are  entered,  press  either  EXECUTE  or  STOP  to  continue  program 
execution.  The  case  identifier  will  be  printed. 

b.  The  display  will  ask  for  the  gain.  This  gain  will  be  used  as  the  forward  loop  gain;  so 
caution  must  be  used  when  interpreting  the  results  of  an  isolated  loop  equivalent  (see  Section 
III). 


c.  The  numerator  and  denominator  coefficients  of  the  closed-loop  transfer  function  are  printed 
out.  Following  this  step,  the  calculator  will  appear  to  be  inactive  for  a long  period  of  time  while 
the  inverse  Laplace  transform  is  calculated. 

d.  The  inverse  Laplace  transformation  will  be  printed. 

e.  The  display  will  request  the  MAX.  TIME  ON  PLOT.  The  time  axis  of  the  plot  is  six 
inches  long,  so  if  the  user  wants  the  tic  marks  on  the  plot  to  coincide  with  the  divisions  on 
graph  paper,  the  maximum  time  specified  should  be  some  multiple  of  six.  After  the  maximum 
time  is  entered,  the  calculator  will  appear  to  be  inactive  for  a long  period  of  time  while  the 
system  response  is  calculated  for  the  time  interval  which  was  specified. 

f.  After  the  maximum  time  is  specified,  the  user  should  insure  that  the  plotter  is  set  up. 

g.  The  printer  will  print  out  some  information,  and  the  response  will  be  plotted. 

h.  The  display  will  then  ask  if  a new  gain  is  desired,  a new  time  range  is  desired,  or  neither 
of  the  above.  If  the  response  to  the  question  is  "N",  the  control  of  the  program  will  be  returned 
to  the  SAP  executive  program.  If  a new  gain  is  requested,  i.e.,  "G",  the  program  will  return 

to  step  b.  If  a new  time  scale  is  specified,  i.e.,  "T",  the  program  will  return  to  step  e. 


^The  term  "steady-state  value"  means  that  the  term  bg/aQ  is  defined.  The  user  should  note 
that  when  the  system  is  unstable,  there  is  no  steady-state  value  even  though  the  term  bg/aQ 
is  defined.  The  program  simply  bases  its  decision  on  the  presence  of  a steady-state  value  on  the 
existence  of  a defined  b^/aQ  and  not  on  the  stability  of  the  system. 
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3.  EXAMPLE  PROBLEMS 


a.  Example  1 

The  time  response  for  the  system  represented  by  the  block  diagram  in  Figure  4 was 
run  for  several  representative  gains.  Figures  22  and  23  are  the  output  and  plot,  respectively, 
for  a forward  loop  gain  of  one,  which  according  to  the  root  locus  will  be  a very  lightly  damped, 
low  frequency  oscillation.  Figures  24  and  25  are  the  output  and  plot  for  a gain  of  6.5.  Figures 
26  and  27  are  the  output  and  plot  for  a gain  of  20,  which  causes  the  system  to  be  unstable. 

The  system  does  have  a steady-state  value  of  the  response,  so  the  responses  are  expressed  as  a 
percentage  of  the  steady  state  value.  Note  that  the  scale  factor  on  Figure  25  is  different  than 
the  scale  factor  on  Figures  23  and  27. 

b.  Example  2 

The  time  response  for  the  transfer  function 


S_±_,1 

$2  (S  + 2)  (S  + 3) 

was  run  to  illustrate  the  format  of  the  plot  for  the  case  where  a system  has  no  steady-state 
value  for  the  time  response.  Figure  28  is  the  output  of  the  TR  subroutine  for  this  example, 
and  Figure  29  is  the  plot.  The  response  is  expressed  in  units,  and  the  scale  factor  is  such  that 
each  tic  mark  on  the  y-axis  is  one  unit.  The  maximum  response  occurs  at  six  seconds,  and  the 
value  of  the  maximum  response  is  3.11  units  as  listed  on  the  output  (Figure  28)  and  shown  on 
the  plot  (Figure  29). 

4.  FLOW  CHART 

The  flow  chart  for  the  Time  Response  subroutine  is  shown  in  Figure  30.  The  complete 
program  listing  for  the  subroutine  is  contained  in  Appendix  F. 
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Figure  27.  TR  Subroutine  Plot  lExamole  1,  Gain 
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Figure  30.  TR  Flow  Chart 
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APPENDIX  B 


SAP-TI  PROGRAM  LISTING  | 
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it  45  60  TO  l?0Oi 

1 700  V-  (SC  1 ? 1 ;.iTTt<Z8“i  > .in  Z8  l'*04  > 

1705  FuR  [“J  TO  C8 

1710  lF  HGSri : 6»  i i+T'-  226  niEtl  1,25 

i715  rii=E: IP unis.  J ]4T)*22t..> 

1 720  Ll'TO  IZ  'iO 
i 7cJ5  Ii  I =EXP  ',  ! L 3 ? I ] i'-T  > 

! 730  1:"  T1  1 , I ili  3 THEN  1745 

i 735  i =T  t ■■  1 L 5 ? I 3-COS  H [ 4 1 j '1*1 1 1 ri  b j I J*S  I N at  4 ^ 1 ] •.  I L I 
1 740  CO  TO  17:.n 

1 v=Yf  i:  Tl  1 ] TI  6>  I J*TfTl  7j  I J-T-TT^Dl 

1750  liEXf  1 

1755 

1760  nC1^  J]-T 
1765  01  2-  !]=■:' 

’ 7701  IF  HH.'  , . [ 4 .1  ' ] ' V]  I HEN  1 7y0 
::V5  ':'i=HBS' 

1780  NEXT  T 
1785  PRINT 

1780  PR' INF  -fiPlX  'nriE  ON  Pl.OT^'  l] 

1795  -PINT  '"'P'^OLUl  E VRLJJE  OT  1 I1hX=-Y1 
1800  .r  SI  2^  i l^-O  TFR  N 1825 
1805  Y0  =SL  1 H : j.  SI  2<  1 ,1 
1810  R7=l 

1815  PHIi  i'i  " -TEHi  .'  CO  IE  YHl.JE  OF  Y='Y0 
1820  CijTO  1040 

1325  RHInT  '■no  STEhIjY  SThIE  VNLUL-- " Z'S- 1 FREE  INTEGROfOR 
r-30  '0*-l 
1835  R '''  = 2 
1840  PPINi 

1345  IF  ,0110  .HEli  185'? 

1850  Y0=1 
1855  ','0=088 (.to:’ 

.860  FOP  M=0  '|-i  Y 
1865  .18=0, 5b1  ONI 
1 8 7 ti  R 0 R 1 = 1 1 1.1  '1 
^'375  = . i 
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1:;:00  IF  VI  1. . I-'W  THEN  1 "'HV 

1035  he, VI  I 

1390  HE  Vi  H 

1395  Y0=3-':V1  '3 

1900  GOTO  19-0 

1905  Y0-J9*Y0 

1910  GOSUG  3465 

1915  FOR'  [ = 1 I'n  61 

1 9 2 0 P L Cl Q L 1.  j 1 ] * 6 0 ■['  1 j C Q [ ? - T ] * ? 0 ' 'i'  i"i 
1925  HEX!  I 
1930  FEH 

1 9 3 5 BIS  P "HE  H G F I H ( G :■  j T I H E S C li  I ..  E < T I' » 1 1 0 C N > ' ? 

1940  INPUT  114 

1945  IF  i:i4  = "G"  THEN  1075 

1950  IF  ]J:f  = "r  THEN  1655 


1955  RETURN 


1 960 

EHB 

1965 

RiEH**** 

•1“  -H"  ■■I- 

■h  -H-  •+■  '•!“ 

•.•r 

•+.  .j:.  h|;  .4..  ;i;,  ..i. 

1970 

REM* 

POLYNOMIRL  ROOT 

3LIBR0IJTINI 

1975 

REM**** 

+• 

4.  4:  .4;. 

+’  +•  •+  rl”  "K  1-  • 

1900 

R1=0 

1985 

R=0 

1 990 

= 0 

1995 

IF  U<2 

THEN 

2435 

2000 

IF  S[3i 

1 

]#3 

THEN  20 

50 

2005 

X ~0 

2010 

Y = 0 

2015 

T=.:i 

2020 

GO 3 UP  2 

G 

70 

2025 

0 " - 0 - 1 

2030 

FOR  I=i 

1 0 0 

2035 

SL3J  I 

s 

r 3)  I 

+ 1 ] 

2040 

HEXl  1 

2045 

GOTO 

CTj 

2050 

If-  3#0 

THEN 

2065 

2055 

R = 0 

2060 

GOTO  20 

i' 

2065 

R=R/S 

2070 

S-- 1 /S 

2 0 7 5 

IF  Pi=a 

THEI-i 

2 8 G 3 

2080 

R 1 = 1 /■R  1 

2085 

S 1-5-*  1 0 

i- 

( - 1 0 

) 

2090 

S[  4 ii  0+  1 

I 

=S[  3 

j 0+ 1 j 

2095 

SL  5?  0+1 

1 

.1 

= 3[  3 

n 0+  1 ] 

2100 

Sr  4j  0 + 2 

J 

* 0 

2105 

3 [ 5 f U + 2 

..i 

= 0 

2110 

FUR  12= 

1 

TO 

20 

2115 

FOR  11= 

0 

TO 

1 STEP 

■1 

2120 

sr  4.  in 

3[  3 ^ 

n ]+pi* 

3f:  4 . 

1 1 + 1 ■] 

4 tr 

i..  1 -...I 

r-  r i-r  T 

T 

9[.i. 

in-:  01  ■■ 

. • - I"  *T*’ 

•J  " 

r j + •'  ■’ 

i 

I 


£130  NE;:'!"  II 

2135  IF  '•  R & S S [ 4 ? 1 ] . • S t 3,11 ) - S 1 > > 0 T R E l-l  2 

2.140  Y=0 

2145  Y^RJ 

2150  T=1 

2155  GOSUb  2.;, 70 

2160  0=0-1 

2165  FOR  1=1  TO  0+1 

2170  313^  I ]=3[4,  I + l ] 

2175  NEXT  I 

2180  GOTO  2050 

2135  IF  315j2]#0  THEN  2200 

2190  R1=R1+1 

2195  GOTO  2205 

2200  Rl=Rl-3[  4"  1 5?  2:] 

2205  FOR  1=0  TO  1 STEP  -1 
2210  S[4i.  I :]=3[3j  I T-R  + SC4H  I + l T--S+3C4J  M2:] 
2215  SCSs  I 1=3[4»  :i  l-R  + 315)  I + l J-S  + S[5-  1+21 
2220  NEXT  I 

2225  IF  S[::3^2  3#0  THEN  2240 

2230  IF  ■:hBS(S[4?2]/SC3M  i:>-31.:'.:-0  THEN  2 

2235  GOTO  2245 

2240  IF  ■ RBS’':SL4j2]/S[3f2j::'-Si::i:'>@  THEN  2 
2245  IF  (flES(S[  4?  1 ]/SC  3j  1 M-Sl <=  0 THE 
2 2 5 £i  T = 3 1 5 ? 2 ..i  ■ ':i  L 4 j 2 1 
2255  1.1=31  5.3 't?-! +SC  5 ? 4 :] 

2260  IF  O#0  ’HEN  2280 
2265  R=R-2 
2270  3=3+ ^ 3+1.:' 

2275  GOTO  2290 

2230  R = F'  i ■ . ;3  [ ^4 . 2 J ■+  S [ 5 » ::3  J - S [ 4 n 13+315.4  3 " U 
2235  3=3  + ''-3!;  4, 2 l + T + SI  4.  1 3 + 31  5.3  3::'/U 
2290  NEXT  12 
2295  3 1=3 1 -:- 10 

2300  IF  3 1 #5U0-M3-6::'  THEN  2110 

2305  PR  I N T " T OL  EF  riNCE=5+ 1 0 ( --6 ') . NO  ROOTS " 

23  1 0 3 i 0+' 

2 :3 1 5 G = R 4 2 4 * 3 

2320  IF  G<0  +NEN  2:365 

2 :3  2 5 X = - P • ' 2 + 3 Gl  R < G 2 

2330  V=0 

2335  T=1 

2340  GOSUB  2670 

2345  7 = - F;  •'  2 - '3  0 F i'  G . ■'  2 

2350  r=i 

2355  GOSLIB  2i:'’’0 

2360  GOTO  2ji35 

2365  X=-R  -2 

2370  V=SQR'  -G::' - ^ 


2390 

2395 

2400 

2405 

2410 

2415 

2420 

2425 

2430 

2435 

2440 

2445 

2450 

2455 

2460 

2465 

2470 

2475 

2430 

2485 

2490 

2495 

2500 

2505 

2510 

2515 

2520 

2525 

2530 

2535 

2540 

2545 

2550 

2555 

2560 

2565 

2570 

2575 

2580 

2585 

2590 

2595 

2600 

2605 

2610 

2615 

2620 

2625 


GOSUB  2670 
0=0--2 

IF  0=0  then  2455'  • 

FOR  I=Oil  TO  1 STEP  ~1 
SL3>1  ]=S[4)  I+2:j 
NEXT  I 

IF  0>2  THEN  2050 
IF  0<2  THEM  2435- 
R=SI:3jO]/S[3jO+1  ] 

S=SC3jO”1  ]/SC3jO+n 
GOTO  2315 

X=-S[ 3? 0 ]/S[ 3? 0 + f ] 

V=0 

T=1 

GOSUB  2670 

RETURN 

END 

REM****  + TIME  RESPONSE  RXES  SUBROUTINE**#**- 

R E f'1  *'******'****+**************************  * 

SCALE  --20  j 80  j -35  j 33 

XRXIS  0ii10ii0j60 

Y R X I S 0 j 5 j - 3 0 j 3 0 

LRFJEL  ( * j 1.5?  1 , 7 ? 0 ) 7 / 1 01 ) 

FOR  1=2  TO  6 STEP  2 
PLOT  10*1 ?0?1 
CPLOT  -3. 6? -1.3 
LABEL  (2525>a/6>*Tl 
NEXT  I 

FORMRT  1F7.2 

FOR  I=-l  TO  1 

PLOT  0<I*20?1 

CPLOT  -9. 8? -0.3 

IF  F7=2  THEN  2565 

LABEL  (2555)1*09*100 

FORMRT  IF 9.0 

GOTO  2575 

LABEL  (2570) 1*Y0 

FORMRT  IE 9, 2 

NEXT  I 

L R B E L (.  •*  ? 2 ? 1 . 7 ? 0 ? 7 . 10) 

PLOT  30?"30?1 
U P L U 1 - 4 . 8 ? - 0 . 6 
LABEL  (*)"TIME  (SEC)" 

PLOT  30?35?1 
B=LEH(Tf) 

H=-<  ■::b-h5)--'2)-*0.  2 
CPLOT  R?-0„6 
LABEL  (*)"CRSE  "Tt 
PLOT  -19?  8?  1 


APPENDIX  G 


INVERSE  LAPLACE  TRANSFORM 
THEORY 

This  appendix  explains  the  theory  that  was  used  to  generate  the  inverse  Laplace  transform 
in  the  Time  Response  Subroutine. 

The  given  system  is  described  by  the  transfer  function, 


Y(s) 


bmS 

_ m 


m 


*■■■■*'^’’*‘’0  (n>m| 


R(s)  apS"  +....  + a^s  + a^. 


(G-1) 


or 


Y(s)  [aps"  + ....  + a-jS  + ap]  = R(s)  {bp^s’’^  + . . . . + b^s  + bp] 


(G-2) 


Taking  the  inverse  Laplace  transform  of  both  sides  yields, 


apy^^^lt)  + + aiy<l>  (t)  + apy(t)  = bp,r<'">(t)  + + bir<1)(t)  + bpr(t)  (G-3) 


where 


Alt) 

dt*^ 


Although  Equation  (G-3)  is  valid  (by  definition)  for  all  inputs  and  all  initial  conditions,  certain 
problems  arise  when  we  attempt  to  take  the  derivative  of  a step  input  at  t = 0.  Thus,  if 
r(t)  A u(t)  and  the  bj's  (i  > 1)  are  not  all  zero,  then  Equation  (G-3)  needs  additional  inter- 
pretation and  clarification.  This  particular  case  will  be  discussed  in  detail  later,  along  with  a 
derivation  of  the  solution. 


The  general  solution  to  any  order  linear  differential  equation  with  constant  coefficients 
is  given  by 


y(t)  = y^lt)  + yp(t)  (G-4) 

\where 

yh(t)  is  the  solution  to  the  homogeneous  equation,  i.e.,  apy^"^(t)  + . . . . -k  apy(t)  = 0 
and 

yp(t)  is  a unique,  particular  solution  which  is  linearly  independent  of  y^lt)  and  is 
non-zero  if  the  right  hand  side  of  equation  (G-3)  is  non-zero. 
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Consider  first  the  general  solution  to  the  homogeneous  equation.  The  general  solution  for 
single  or  multiply  repeated  real  roots  can  be  expressed  as 


I J 

= Z)  S Cjj  tj'^  e^jt 

i=1  j=1 


(G-5) 


where 


Aj  is  the  distinct  root  of  the  characteristic  equation, 

I is  the  total  number  of  different  roots^ 

J is  the  multiplicity  of  the  i^^  root,  i 
Cjj  are  the  undetermined  coefficients. 

The  general  solution  for  complex  roots  can  be  expressed  as 


I 


~ S cos  wjt  + B,j  sin  Wjt]  (G-6) 

i*1  / 


where 

/ 

a i is  the  real  part  of  the  i^^  cr/mplex  root  pair, 
wj  is  the  imaginary  part  of  the  i^^  complex  root  pair, 

1 is  the  total  number  of  cbmplex  root  paris, 

Aj,  Bj  are  the  undetermined  coefficients. 

Therefore,  the  general  homogeneous  solution  is  the  sum  of  Equations  (G-5)  and  (G-6).  The 
case  of  multiply  repeated  coqiplex  root  pairs  is  so  rare  and  cumbersome  that  it  is  omitted 
from  the  computer  program. 

Next  consider  the  particular  solution  for  the  differential  equation.  The  driving  function 
tor  the  Time  Response  Subroutine  is  a unit  step  applied  at  t = 0.  The  particular  solution  must 
satisfy  the  equation 


^nyp'"^<^)  + . . . . ■^  a^ypit)  = b^uit)  + b^Sft)  + ....-*•  bf„  6 ‘ '*ht) 


(G-7) 


where 


u(t)  is  the  unit  step  function  and  S^'^^ft)  is  the  q^^  derivative  of  the  delta  function. 


which  is  the  (q  + 1)^^  derivative  of  the  unit  step  function. 

The  problem  with  the  unit  step  is  that  it  is  undefined  at  t = 0.  Evaluating  Equation  (G-/ ) at 
t = 0'*’  will  alleviate  the  problem  of  an  undefined  unit  step,  and  also,  the  delta  function  and 
its  derivatives  will  be  zero.  However,  evaluating  the  equation  at  t = 0^  means  that  we  will 
have  to  calculate  new  initial  conditions  for  the  system  since  the  system  was  subjected  to  an 
infinite  excitation  at  t = 0 and  thus  the  system  at  t = 0^  will  not  be  at  rest. 

Therefore,  Equation  (G-7)  becomes 

anyp(n)(0+)  + ....  + a^yplO+l  = b^  (G-8) 


where 

u(0+)  = 1 

6<q>(0‘^)  = 0 for  q >0 
Assume  a particular  solution 

yp(t)  = (G-9) 

where 

q = the  subscript  of  the  coefficient  of  the  lowest  order,  non-zero  derivative  in  the 
left-hand  side  of  Equation  (G-8). 

This  general  form  of  the  particular  solution  is  used  because  it  is  also  valid  for  systems  which 
have  one  or  more  poles  at  the  origin. 

Differentiation  of  Equation  (G-9)  yields 


(G-10) 
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yp(t)  = Kt^ 
yp<l>(t)  = qKt^'' 

• 

= q!  K 

yp<‘’  ^>(t)  = 0 

yp<">(t)  = 0 


A 


VI 

t 


' i 


r 


When  the  set  of  equations  in  (G-10)  is  substituted  into  Equation  (G-8)  and  the  resulting 
equation  is  evaluated  at  t = 0"^  we  have 

a qlk  = b 
q o 


or 


Therefore, 


yp(t)  = 


q! 


(G-11) 


(G-12) 


Substituting  Equations  (G-5),  (G-6),  and  (G-12)  into  Equation  |G-4),  differentiating  the  result 
(n  - 1)  times,  and  evaluating  at  t = 0"^  yields  a set  of  equations. 


But 


and 


y(0+)  = yp(0+)  yhr(0+)  -i-  yhc(0+) 

y(n  - 1)(o+)  = yp(n  - 1)(o+)  + y^^(n  - 1)(o+,  + yjn  - 1)(o+, 
Vn^'^No'*')  = 0 for  k ifc  q 


yp^'^^(0''')  = °/^q  for  k = q 

Therefore, 

yhr<o‘'»  + yhc<o‘')  = y<o'') 

yhr<^><0+)  yhc^^’^iO'^)  = y<‘^>(0+)  - 


yhr‘"  ‘ ^^0+)  -t-  yhc<"  ■ ’>(0+)  = y<"  ' 1){0+) 


(G-13) 


(G-14) 
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where,  in  general,  the  y are  the  initial  conditions  of  the  system  at  t = 0''^. 

The  final  step  prior  to  solving  the  set  of  equations  (Equation  G-14)  for  the  undetermined 
coefficients  associated  with  the  homogeneous  solutions  is  to  determine  the  initial  conditions 
at  t = O'*'.  Since  we  know  that  the  initial  conditions  at  t = O'  are  all  zero,  we  should  be  able 
to  compute  the  state  of  the  system  at  t = O'*'.  Rewriting  Equation  (G-3)  with  r(t)  = u(t),  yields 


a^y^'^^ft)  + -r  agylt)  = b^uft)  + bi6(t)  + + br„6<'"  ' Dft) 


{G-15) 


Let's  integrate  this  equation  n-times  from  t = 0'  to  t.  The  first  integration  yields 


'n  /v 


an  y y ^'^'(t)dt  + ....  + a^J  y(t)dt  = I u(t)dt  + b^  /5(t)dt  + . . . . 
O'  0-  •'o-  *^0- 


,f  y(t)dt  = ^Q  f u(t)dt  + b■^  fiix 

n- 


(G-16) 


+ br„  - l){t)dt 


Integrating  Equation  (G-15)  n times  will  yield  n equations  for  the  n derivatives  y(t), 
y(t),  . . . y"  _ '(t).  Evaluating  these  n equations  a t = 0"^  will  result  in  the  n initial  condi- 

tions yfO'*'),  y(0‘‘'),  . . . y"  ■ ' (O"*^).  For  example,  the  equation  which  results  from 
integrating  Equation  (G-15)  k times,  when  evaluated  a t = O'*"  yields 

a^y(n  - k)(Q-fj  + + g^y  (-k)(o+)  = boU<''^)(0'^)  -i-  5^6  <-k)(0+) 

(G-17) 

. . . . -H  bgS  ■ '^)(0'^) 

where,  in  general, 

f(q)(0+)A  f(t)  I t = o^  forq  = 0 
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f<q)(o+)  A.  f(t)  . 

' t = 0^ 


for  q > + 1 


0"^  t t 

f(q)(0+)  4 Iff  f|t)dt  = ) for  q < -1 

0 "O'  O' 

q times 
and 

f(q)(0  ) = 0 due  to  the  fact  that  the  system  was  at  rest  at  t = O'. 

Several  observations  can  be  made  about  the  above  equation. 

(a)  Since  y(t)  = 0 at  t = O'  and  must  be  finite  at  t = O’*",  then  all  integrals  of  y(t)  evaluated 
from  t = O'  to  t = O’*”  must  be  zero,  i.e., 

Y^q)(0+)  = 0 for  q < 0 

(b)  The  delta  function  terms  on  the  right-hand  side  are 


6<^>(0+)  = 1 

for  q = -1 

6<q>(0‘^)  = 0 

for  q >0 

5(q)(0+)  = 

1 =0  for  q<  -2 

lt  = o-^ 

(c)  The  unit  step  function  on  the  right-hand  side  is 


ulQ-")  = 1 


for  q = 0 
for  q <-1 


Since  the  derivatives  of  the  step  function  are  expressed  in  terms  of  the  & - function,  the  case 
for  q >1  is  not  considered  since  no  derivatives  of  u(t)  will  appear  in  any  of  the  (n  ^ 1)  equations. 


(d)  Since  evaluating  (q  > 1)  involves  a knowledge  of  ' ^^0'*’),  y^*^  ' ^*(0'*’) 

. . . . y(O^),  the  proper  procedure  for  evaluating  the  n-initial  conditions  involves  first  evaluating 
ylO'*')  using  Equation  (G-17)  for  k = n.  Then  y(0‘*’)  can  be  evaluated  using  (G-17)  with  k = n 1 
and  the  previously  computed  value  of  ylO"^).  This  procedure  is  continued  until  y’^'MO'*')  is 
computed  from  (G-17),  where  k = 1. 

Considering  these  observations,  evaluating  Equation  (G-17)  for  several  values  of  k yields 
y(0+)  = 0 


y(n  - m - 1)|g-l-|  = 0 
y(n  - m),o-t-)  = 


y(n  ' m + 1)(o+)  = — !—  b^  . T - a^  . y " ' '^(0+) 

^n  [ J 

y(n  - 1)(o+)  = __L  bi  - a^  . 1 y<"  ■ 2)(0+)  ....  -aiy(0+) 

^n  1_ 

The  set  of  equations  (G-18)  can  be  written  in  a general  form  as: 


(G-18) 


,ln  - m + ql|o+|  . ^ ^ ^ 

" L k^ 

for 

0 ^ q m - 1 

As  shown  in  the  set  of  equations,  the  value  of  the  initial  conditions  for  any  derivative  of  order 
lower  than  (n  - m)  is  zero,  i.e., 

y(n  - m -H  q)(o+)  = o for  q < 0. 
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an  . kv'"  ■ ^ ■ ''*(0^)  (G-19) 


Once  the  initial  conditions  have  been  evaluated  using  the  above  recursive  relationship,  the 
final  step  involves  determining  the  remaining  constants  (coefficients)  when  the  set  of  equations 
(G-14)  is  evaluated  at  t = 0'*'.  The  most  convenient  method  to  accomplish  this  is  to  use  the 
following  matrix  representation  for  n-simultaneous  linear  algebriac  equations  with  n unknown 
coefficients: 


‘All  • • ■ • 

A-] 

'n 

■cr 

-y(0+)-yp(0+) 

• 

• 

• 

• 

• 

• 

• 

_ 

• 

• 

• 

• 

• 

• 

• 

• 

• 

_^n1 

Ann. 

■ ■'>(0+)-yp<"  - 

(G-20) 


or 


[A]  [C]  = [b] 


where  the  column  matrix  [C]  contains  the  undetermined  coefficients.  The  coefficient  matrix 
can  be  solved  as  follows: 


[C]  = [Aj-1  [b] 


After  all  of  the  coefficients  have  been  determined,  Equations  (G4),  (G-5),  (G-6),  and  (G-12) 
give  the  general  solution  for  y(t). 
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